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Josephson junctions fabricated on the surface of three-dimensional topological insulators (TI) show a few 
unusual properties distinct from conventional Josephson junctions. In these devices, the Josephson coupling 
and the supercurrent are mediated by helical metal, the two-dimensional surface of the TL A line junction of 
this kind is known to support Andreev bound states at zero energy for phase bias 7r, and consequently the 
so-called fractional ac Josephson effect. Motivated by recent experiments on Tl-based Josephson junctions, 
here we describe a convenient algorithm to compute the bound state spectrum and the current-phase relation for 
junctions with finite length and width. We present analytical results for the bound state spectrum, and discuss the 
dependence of the current-phase relation on the length and width of the junction, the chemical potential of the 
helical metal, and temperature. A thorough understanding of the current-phase relation may help in designing 
topological superconducting qubits and manipulating Majorana fermions. 



I. INTRODUCTION 

Josephson junctions have been the key elements in super- 
conducting devices such as SQUIDs 1 . In the past decades, 
they have also become the staple components for supercon- 
ducting qubits 2 in the general architecture of circuit quantum 
electrodynamics 3 . Recently it was pointed out that Josephson 
junctions patterned on the surface of topological insulators 
(TI) can be used to create and manipulate Majorana fermions 
for topologically protected quantum computation 4 . More gen- 
erally, Tl-based topological qubits can be integrated with the 
standard superconducting qubits to achieve a new hybrid plat- 
form for information processing 5 6 . These prospects motivate 
us to carry out a detailed investigation of the equilibrium prop- 
erties of Tl-based Josephson junctions. 

A fundamental property of any Josephson junction is its 
current-phase relation (CPR), /(</>), where / is the equilib- 
rium supercurrent through the junction and <fi is the supercon- 
ducting phase difference across the junction 7 . For the well 
known tunnel junction originally considered by Josephson, 
the current-phase relation is simply I(<j)) = I c sin cj>, with I c 
being the critical currenP. By contrast, the CPR of a pin- 
hole junction, also commonly referred to as a superconducting 
constriction, has a rather different form, I(<j)) = I c sin((/>/2) 9 . 
The CPR for junctions between unconventional, such as p- 
wave or <i-wave, superconductors are generally more compli- 
cated (for review, see Ref. |7]). The CPR is sensitive to the 
pairing symmetry of the superconductor as well as the micro- 
scopic scattering channels and amplitudes, and can be mea- 
sured directly in experiments. Anomalies in the CPR often 
point to new physics. The CPR also controls the dynamic 
properties of the junction, especially when the phase dynam- 
ics is slow compared with the inverse gap. Understanding the 
CPR is thus important for designing superconducting circuits 
and qubits. 

The main problem we address is how to find the CPR for 
the Josephson effect mediated by a new state of matter, the 
helical metal at the surface of three dimensional topological 
insulators Helical metal, consisting of massless Dirac 



electrons with spin-momentum locking, is much more exotic 
than graphene. It is only "a quarter of graphene," with an odd 
number of Dirac cones (for simplicity we only consider a sin - 
gle Dirac cone at k = as found in Bi2Se3 and Bi2Tes!^^). 
Microscopically, the supercurrent flow is tied to the process of 
Andreev reflection, as in the well known case of the Joseph- 
son effect through a two-dimensional electron gas. However, 
the Andreev reflection of helical Dirac electrons differs from 
that of conventional electrons with quadratic dispersion. One 
then expects that new scattering kinematics such as specular 
Andreev reflection, which was discovered in the context of 
graphene by Beenakkei 12 , will strongly influence the Andreev 
bound state spectrum and consequently the CPR in certain 
regimes. For example, we will present an interesting scaling 
relation between the critical current and the length of the junc- 
tion which is unique to helical metal with chemical potential 
right at the Dirac point. In this case, the supercurrent may be 
thought as being carried by evanescent waves, but it does not 
decay exponentially with the length of the junction. 

Many of the new features of the Josephson effect through 
helical metal were recognized and discussed in the pioneer 
work of Fu and Kan&H Most notably, they discovered that 
a short line junction at phase bias tt features a linearly dis- 
persing Andreev bound state spectrum with a robust crossing 
at zero energy for transverse momentum k y = 0, so the line 
junction is a "Majorana quantum wire' 4 . The objective of this 
paper is to generalize their analysis to junctions with finite 
length and width, and systematically investigate the effects 
of finite chemical potential of the helical metal and tempera- 
ture. The motivation is to make predictions that can directly 
compare with experiments. Finding the CPR for such junc- 
tions turns out to be algebraically cumbersome. We outline a 
procedure that is conceptually simple while straightforward to 
implement. This also enables us to find a few new analytical 
results for finite size junctions. 

Several groups have successfully fabricated Josephson 
junctions of various lengths on exfoliated flakes or epitaxial 
thin films of Bi2Se3 and observed supercurrent^. It re- 
mains unclear that the supercurrent is entirely due to the TI 
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surface states in all these published results, because the TI 
bulk also conducts for many samples used in experiments. 
This problem can, however, b e circ umvented by applying a 
back gate 19 21 , chemical doping 19 22 , or adopting the new gen- 
eration of so-called ideal topological insulators 23 25 where the 
chemical potential is tuned inside the bulk gap. Thus, we will 
focus on the physics associated with the helical metal, and as- 
sume conduction through the bulk has been eliminated using 
one such technique. 

The current-phase relation for Josephson junctions on the 
TI surface has been investigated theoretically with ferromag- 
nets sandwiched between the superconductor^ 2 ^, which in- 
troduces an energy gap to the helical metal. In Ref. 120 the 
Josephson effect through helical metal was considered, but 
the superconductors are assumed to be conventional BCS su- 
perconductors, which differ substantially from the Fu-Kane 
model 4 adopted here. The anomalous Josephson current via a 
vortex pinned to a hole drilled through a TI slab was studied 
in Ref.|29] 



II. MODEL HAMILTONIAN AND SOLUTION STRATEGY 

The Josephson junction under consideration is shown 
schematically in Fig. [T] The chemical potential of the topo- 
logical insulator is assumed to be tuned inside the gap so there 
is no bulk conduction. The two-dimensional surface of TI, the 
helical metal, is modeled by the HamiltoniarP^^ 



/i M (k) = v(a x k y - (j y k x ) - hm- 



(1) 



Here, a XyV are the Pauli matrices in spin space, k = (k x , k y ) is 
the two-dimensional surface momentum is set to 1 through- 
out the paper), and v is the velocity of the helical Dirac elec- 
trons. As argued in Ref. 4, the presence of an s-wave super- 
conductor (S) induces a pairing interaction between the helical 
Dirac fermions at the surface of the topological insulator, and 
gaps out the surface spectrum. The S-TI interface can then be 
modeled elegantly by a simple matrix Hamiltonian in Nambu 
space (we follow the convention of Ref.fTTT). 



h s (k) 
-icr 7/ A* 



ia y A s 



where 



h s (k) = v(cr x ky - (j y k x ) - us- 



(3) 



In general, we allow the chemical potential /is and /im to be 
different. For example, one can add gate control over hm in 
the helical metal region. The model Eq. ^ has been shown to 
be accurate at low energies for both weak and strong coupling 
between S and TI using self-consistent calculations^"^. 

The whole system is translationally invariant in the y direc- 
tion if the width of the junction is infinite, W — » oo. We define 
x = to be on the boundary between the left superconductor 
Si and the helical metal, and x = L to be the boundary be- 
tween the helical metal and the right superconductor S2. The 
phase of Si is chosen to be 0, while the phase of S 2 is denoted 




FIG. 1: (Color online) a) Schematic of a Josephson junction through 
helical metal, b) The particle (e) and hole (h) branch of the excitation 
spectrum for the helical metal, c) The gapped spectrum (solid lines) 
of Hs describing the superconductors. The dashed lines show the 
particle and hole excitations if the proximity effect is absent (A = 0). 
\im and fis are measured from the Dirac point, and are different in 
general. 



<j). The Hamiltonian is piece- wise constant in each region and 
has the following generic form: 



«(k) = 




(4) 



where k± = k x ± ik y , fi = hm for < x < L and /is 
elsewhere, A(x > L) = A e i(f) , A(x < 0) = A , and A(0 < 
x < L) = 0. A is the bulk gap of both superconductors. 
Note that retaining the hole sector in the helical metal region 
is crucial for our discussion of the Josephson effect. 

To find the CPR, we solve for the Andreev bound state 
(ABS) spectrum, i.e., solutions to the eigenvalue problem 



H^ = E^ 



(5) 



with \E\ < Ao for given phase difference and k y . While the 
problem is conceptually equivalent to finding the bound states 
in a square potential well, the algebra is more complicated 
because of the matrix structure of H. Our basic strategy is to 
search for those values of E for which the eigenvector i/j is 
continuous at x = and L. 

To this end, it is important to known all the traveling as well 
as evanescent wave solutions of Eq. ^ in each region. We 
(2) use a triclP^ to organize these solutions which turns out to be 
crucial in obtaining the ABS spectrum for general parameters. 
We first rewrite 1-L as 



/0100 
10 
1 

Vo 1 



vk x 



I 



-fi 
-ivky 


A* 



ivk, 



y 
-/i 

A* 





-A 



ivk v 



A 


-ivky 



and then rearrange Eq. ^ into an eigenvalue problem for k x , 

}C?p = k x ijj, (6) 



where the matrix JC is non-Hermitian, 
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(7) 
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and the eigenvalue k x is generally complex. For given cj), E, 
and k y , we can solve Eq. ^ easily. The complex k x solutions 
describe evanescent waves. Physically, the wavefunctions for 
ABS decay inside the superconductor. As another example, 
the bound states for /im = involve evanescent rather than 
traveling wave solutions of the helical metal Hamiltonian Km 
(and its counterpart in the hole sector). 

To demonstrate the procedure, we start by discussing the 
simple case of k y = analytically. This corresponds to 
Josephson coupling through a TI nanoribbon with small W, 
where transverse quantization of k y only allows the single 
mode k y = below the energy scale Ao. Note that the length 
of the nanoribbon L is kept general. 



parts of k x for S2, and similarly the two solutions with nega- 
tive imaginary parts for Si. All four eigenfunctions are used 
for the helical metal. Requiring continuity of ip at x = and 
x = L, we can solve for the eight unknown coefficients. Do- 
ing so, it is straightforward to show that the unique solution 
for zero energy bound states is always at 

(j) = 7r. (13) 

This agrees with the results of Ref . |4] 

B. Finite Energy Solutions 



A. Zero Energy Solution 

We first examine the case of E = and k y = 0. Solving 
([6]), we find eigenvalues of 



k\ = -(iAo+fJLs) 

k 3 x = ^(-iAo + fjbs) 
k x = ^(-iAo-fis) 
with associated eigenvectors (not normalized) 

fa = (ze^,ie^,-l,lf 
^ 2 = (-ie^ie^,l,lf 
tl>3 = (-ie i +,-ie i *,-l,l) T 
fa = (ie^,-ze^,l,lf 

for the superconductor, and degenerate solutions of 

k x = k x = ~L 1 m 
1 



with 



fa = (1,1,0, Of 
</> 2 = (0,0,-1, if 

fa = (-1,1,0, of 

fa = (0,0,1, if 



(8a) 
(8b) 
(8c) 
(8d) 



(9a) 
(9b) 
(9c) 
(9d) 



(10a) 
(10b) 



(11a) 
(lib) 
(11c) 
(Hd) 



for the helical metal. The wave function in each region is a 
superposition 



^ = 5>^-e^ . 



(12) 



Now we move on to general bound states at finite energy 
for k y = 0. It is convenient to parameterize E using an angle 
/?, E — A cos /?, for |i£| < A . For the superconductors, we 
find eigenvalues and eigenvectors of 



kl = -{ii s + iA sin/ 

kry. 



-(fis - sin/ 



v 

kl = -(-/us + iA sin/ 

v 

k x = "("Ms - iA sm/ 



and 

^i = (e^+«,e^),-l,lf 
^ = (e < W-«,c < W-«,-l,l) T 
^ 3 = (e^- /3) ,-e^),l,lf 
^ 4 = (e^ +/3) ,-e^),l,lf 

In the helical metal region, we have eigenvalues of 

kl = ^(hm + A cos/3) 
kl = -{hm ~ Aqcos/S) 



v 

kl = -(-Mm - Aocos/3) 
kl = ^(~Mm + Aqcos/S) 



(14a) 
(14b) 
(14c) 
(14d) 



(15a) 
(15b) 
(15c) 
(15d) 



(16a) 
(16b) 
(16c) 
(16d) 



and the same eigenvectors as in Eqs. (JTTJ. 

By matching the wavefunctions at the two boundaries, we 
find that the bound state energy E has to satisfy the following 
transcendental equation, 



£? = ±Aooob[^±|]. 

hv 2 



(17) 



Because the ABS wave function must go to at x = ±00, we 
only include the two eigenfunctions with positive imaginary 



This analytical result demonstrates a remarkable feature of 
these junctions: the ABS energy does not depend on hm or 
lis (this is on ly valid for k y = 0). As a sanity check, for 
E = the solution is 6 = ±7r which we have found earlier: 
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zero energy states are always at (j) = ±tt. In the short junction 
limit, L — »• 0, we have 



£ = ±A cos(|). 



(18) 



which was obtained by Fu and Kane earlier in Ref. 21 For 
long junctions, there are many quantized ABS levels which 
only slowly disperse with (/>, 



nn hv 

TI' 



(19) 



where n is an integer. 



III. ANDREEV BOUND STATES 

Now we describe how the ABS spectrum can be obtained 
numerically for the general k y . Continuity of each wave func- 
tion component at x = and x = L gives in total eight equa- 
tions. These boundary conditions can be organized neatly into 
a matrix equation in the form of Az = 0, where the 8x8 
matrix A is a function of E, (j), and k y , and z is a column 
vector containing the eight unknown coefficients. A nontriv- 
ial solution requires the determinant of matrix A, D = det A, 
to be 0. Then to find the allowed energies E for a given k y 
and (/>, we just have to find the zeroes of D(E) in the range 
— Ao < E < Ao. Given the particle-hole symmetry of the 
problem, we only need to look in the range < E < Aq. 
D(E) is in general complex, so we look for zeroes of its ab- 
solute value. Numerically it is much easier to search for min- 
ima of \D(E) \ rather than the zeroes directly. Therefore, our 
algorithm starts by splitting the energy range < E < A 
into N equal slices. Within each of these slices, we perform 
a standard golden section search for minima, as described in 
chapter 10 of Numerical Recipes 34 . After a point is identified 
to be a local minimum, we check if \D(E)\ at that point is 
close to zero (less than a tiny error tolerance). We also check 
the endpoints of the slices to see if they are zeroes. N has 
to be sufficiently large to exhaust all the zeros. Finally we 
check and eliminate unphysical solutions, e.g., those leading 
to a matrix A of rank lower than 8. 

Fig. [2] compares the ABS spectrum {E n ((f))} for k y = 
and k y = 0.2£ -1 , where the "coherence length" £ = hv/A , 
and fiM = 0. Here we assume the junction is infinite in the y 
direction and k y is a good quantum number. In each case, we 
see that as the junction length L is increased, more branches 
of ABS show up. The key difference is that in the case of 
finite k y , the crossings at <\> = and 7r observed for k x = 
are now replaced by avoided crossings, and E n {(j>) become 
smooth functions. 

We notice two singular features in the function E n ((f)). The 
first is the crossing at zero energy and <p = n. Strictly speak- 
ing, the crossing is only for k y = 0. But for small k y , the 
change in the slope of E n {4>) is still rapid, and even for large 
values of k y , the slope of E n ((j)) changes sign at <p = n. The 
second is the merging of the ABS into the quasiparticle con- 
tinuum at E = ±A at some finite values of <j) which we 
denote <j> c . Both features will lead to a sudden change in the 
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FIG. 2: Andreev bound state energies as functions of the phase dif- 
ference (j). Only E > are shown. Top panel: k y = 0. Bot- 
tom panel: k y — 0.2£ _1 . The junction parameters are us — 2Ao, 
Mm = 0, W = 00, and the junction length L is measured in units of 

£ hv/Ao. 



slope, dE n {4>)/d(j). This, provided that the ABS is occupied, 
will leave fingerprints in the current-phase relation at low tem- 
peratures. 



IV. CURRENT-PHASE RELATION 

After the ABS spectrum {E n (<p)} is found for given k y , 
the k y -resolved supercurrent is given by the phase dispersion 
of{E n }, 



yi 1 



2e ^ dE n 

~h 2-^ rlrh 



1 



(20) 



where T is the temperature, the Boltzmann constant ks is set 
to 1, and the sum is over all ABS energies (continuum quasi- 
particle excitations give zero net contribution to the supercur- 
rent). We can further exploit the particle-hole symmetry to 
rewrite this as 



E n >0 r 



(21) 
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Numerically, we approximate the derivative as 

dE _ E((j) + e)-E((l)-e) 



2e 



(22) 



for small 
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FIG. 3: Momentum (k y ) resolved supercurrent at zero temperature 
as function of (j). Top panel: k y = 0. Bottom panel: k y = 0.2£ _1 . 
The junction parameters are the same as Fig. [2] 



Fig. [3] shows the I(k y , (j)) corresponding to the ABS spec- 
trum shown in Fig. [2] For k y = 0, the sudden jump and sign 
change of I at <j) = 7r can be traced back to the zero energy 
crossings in the ABS. Due to the crossing, the slope of oc- 
cupied ABS at zero temperature experiences a sign change. 
For a short junction, I(k y = 0,0) = J c sin(0/2), in agree- 
ment of the analysis of Ref. 4 The remaining sharp turns, 
i.e., discontinuities in the derivatives of /, occur at (j) c where 
the ABS reaches the superconducting gap and is absorbed into 
the quasiparticle continuum. For finite k y , the sudden drop at 
(j) = 7r is replaced by a smooth variation, but the sign change 
of / remains. Particularly, we observe that I ((f) = 7r) = 0. 

Now we discuss the effect of temperature. Fig. [4] illustrates 
the evolution of I{k y = 0, </>) for a short and a long junction. 
As T rises, the sawtooth-shaped CPR gradually becomes a 
sine function at high temperatures. This is due to the thermal 
population of all ABS levels which tends to smooth out the 
sudden jump at (j) = tt. By comparison we see that the sharp 
turns in I(<j)) for long junctions survive even at finite tempera- 
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FIG. 4: Effect of temperature on I{k y — 0,0). Upper panel: a short 
junction with L — 0.01£. Lower panel: a long junction with L — £. 

lis — 2Ao, \im — 0, W = 00. 



tures. These features are thus in principle observable in future 
experiments. 

In order to get the total supercurrent, we must sum over all 
possible values of k y , 



(23) 



7(0) = / dkyliky^). 



For infinitely wide junctions, the integration goes from — k* F 
to kp, with kp = (fiM + Ao)/v. For a junction with finite 
width W, we assume for simplicity open boundary conditions 
at y = and W (this may not accurately represent certain 
experiments, such as those in Ref. [17]). Then k y is quantized, 
ky = irn/W, with n being an integer, and the integral in 
Eq. (23l is replaced by a discrete sum over \ky\ < k* F . This 
yielasthe current phase relation. Fig. [5] compares the CPR of 
junctions with various widths. As the number of transverse 
modes increases, the total current increases accordingly. The 
overall shape of the CPR remains roughly the same, however. 

Finally, we examine how the chemical potential of the he- 
lical metal, ^m, affects the CPR with us fixed. For k y = 0, 
changing {1m has no effect. This is explained by the lack of 
any hm dependance in the ABS spectrum we found analyt- 
ically for k y = 0, as seen in Eq. (17). For k y ^ 0, how- 
ever, fiM does affect the current. As shown in Fig.|6j the k y - 
resolved current generally increases with hm until hm = Us, 
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FIG. 5: Total supercurrent for junctions with finite width W. /j,s = 
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FIG. 7: Total supercurrent for junctions with various helical metal 
chemical potential \±m- \is — 2A , L — £, T = 0.1A , VK=16£. 
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FIG. 6: Momentum resolved supercurrent for junctions with 
various helical metal chemical potential \±m- Ms = 2Ao, k y — 
0.2£-\L = £,T = 0.1A . 



at which point it saturates and varies very little. The total 
current, as shown in Fig. [7] has a relatively consistent shape, 
but its magnitude increases with j^m primarily because more 
transverse modes are available for large hm- 



0.4 



0.2 



0.2 0.4 0.6 0.8 1 1.2 1.4 

FIG. 8: The supercurrent at zero temperature through a single mode 
(k y = 0) TI nanoribbon of length L. /is = 2A , \im — 0. 



point (for review, see Ref. 35). There it was realized that 
graphene behaves very much like a disordered metal. I c 
reaches a minimum value for \im = 0, where I c scales with 
1/L. For helical metal, we find that the numerical data in Fig. 
[5] can be fit by 



V. SCALING OF THE CRITICAL CURRENT WITH 
JUNCTION LENGTH 

Now we consider the critical current I c for a junction 
through a topological insulator nanoribbon of length L with 
a single transverse mode k y = (due to small W). Fig. [8] 
shows the numerically calculated I c as a function of L. Since 
Mm = 0, there are no states (propagating modes) available at 
the Fermi level of the helical metal. Naively, one would expect 
the critical current to decay rapidly with L — indeed, the de- 
cay should be exponential in L if the helical metal is replaced 
by a normal insulator. This interesting problem was inves- 
tigated for the Josephson effect through another semimetal, 
graphene, with chemical potential tuned right to the Dirac 



(24) 



Actually, since we have derived analytical results for the ABS 
spectrum for k y = 0, Eq. ( [17] ), we can derive the equation 
above analytically. Taking the derivative of Eq. §FJ\ with re- 
spect to <j) on both sides and recognizing that that the critical 
current corresponds to approaching 7r, we obtain Eq. ([24). 



VI. SUMMARY 

We have described how the current-phase relation of the 
Josephson effect through helical metal can be obtained for 
general parameters, including the length and width of the 
junction, the chemical potentials, and temperature, as moti- 
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vated by recent experiments. A few useful analytical results 
that we derived, especially for single mode nanoribbon junc- 
tions, yield nontrivial predictions on the scaling of the crit- 
ical current with the length of the junction. The numerical 
algorithm outlined here can be implemented to model the su- 
percurrent flow in experiments where the current is entirely 
carried by the surface of the helical metal. A detailed under- 
standing of the spectra and CPR of these Josephson junctions 
will contribute to the general goal of using them for supercon- 
ducting devices and topological qubits. 
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